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Ab. stra.c£ 

The MATRICS flow solver calculates the 
inviscid transonic potential flow about a wing/body 
semi -configuration. At present, work is In progress 
to extend MATRICS to take viscous effects into 
account through coupling with a boundary layer 
solver. This solver, MATRICS -V, is based on robust 
calculation methods for the boundary layer, the 
outer wing flow and their interaction. MATRICS -V is 
intended for (Inverse) design purposes. The 
boundary layer and wake are based on an integral 
formulation of the unsteady first order boundary 
layer equations, the inviscid method is the 
existing MATRICS potential flow solver and the 
interaction algorithm is of the quasi-simultaneous 
type. 

The paper gives a progress report on the 
coupled potential -flow boundary -layer method for 
transonic wing/body configurations. 

1. Instruct tan 

Computation methods for three-dimensional 
transonic potential flow are an important component 
In design systems for civil aircraft wing/body 
configurations. Accurate performance prediction 
under these conditions (e.g. lift-drag analysis) 
requires that the transonic potential flow solver 
is coupled with a boundary layer solver to account 
for the viscous effects on the wing pressure 
distribution in a sufficiently accurate way. For 
some years NLR has available its own developed 
system MATRICS (Multi -component Aircraft TRansonic 
Inviscid Computation System) for the calculation of 
the three-dimensional inviscid transonic potential 
flow about wing/body configurations (Ref. 1). The 
MATRICS -derivative MATRICS -V -now under develop- 
ment- will calculate the Interaction of the 
inviscid potential outer flow and a viscous 
boundary layer on the wing of a transonic transport 
ving/body configuration. The ultimate objective of 
the development of MATRICS -V will be the embedding 
of this system in a wing design system, also 
currently under development at NLR. 

This paper gives a progress report on the 
development of MATRICS-V. Firstly, a description 
will be given of the MATRICS three-dimensional 
transonic potential flow solver, being the starting 
point for the development of MATRICS-V. Secondly, 
the requirements for the development of the 
viscous -inviscid interaction solver will be 
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formulated. Robustness of the interaction algorithm 
has been formulated as the main requirement. Next, 
the basic concepts of the viscous- Inviscid 
interaction solver will be given. Subsequently, a 
short system overview will be given of the MATRICS - 
V flow solver, followed by an analysis of the 
boundary layer equation system and on the Inter- 
action law. Finally, some preliminary computational 
results will be presented. 

2. Starting position for the 
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The MATRICS-V viscous - inviscid Interaction 
solver is a follow-up to the MATRICS three- 
dimensional inviscid transonic potential flow 
solver. This flow solver solves the full potential 
equation in strong conservation form on a grid of 
C-H or C-0 topology. This grid is generated using 
the MATGRID grid generator (Tysell and Hedman, 
Ref. 2). The solver uses a fully conservative 
finite volume discretization and a multigrid 
solution method. The discretization scheme is 
second order accurate in the mesh size in subsonic 
parts of the flow, and first order accurate in 
supersonic parts of the flow. For the capture of 
supersonic/subsonic shock waves a Godunov- type 
shock operator is used. Options for fully- 
conservative as well as non- conservative shock- 
capture are available. Details on the MATRICS flow 
solver can be found in references 3, 4. 

MATRICS provides data for a lift -drag- diving-moment 
analysis. Substantial research has been performed 
on the reliable prediction of drag by the MATRICS 
flow solver. Findings of this research have been 
reported in reference 5. 

The development of MATRICS-V is partly based 
on experience at NLR with the modelling of two- 
dimensional strong viscous -inviscid interaction on 
airfoils (Refs. 6, 7). 

1 l Reaulrenents 

MATRICS-V is designed to calculate the 
influence of the wing boundary layer and wake on 
the inviscid potential outer flow about a given 
transport- type wing/body configuration from 
subsonic up to and including transonic cruise 
conditions. Laminar as well as turbulent boundary 
layers should be calculated, where the turbulent 
boundary layer is allowed to be mildly separated. 
The MATRICS-V method should be robust; a guaranteed 
converged solution should be obtained for realistic 
flow conditions . The method should also be about 
ten times faster than a three-dimensional Reynolds - 
averaged Navier- Stokes solver, to enable wing 
design applications. 

4. Basic concents 

The coupling of an Inviscid outer flow with a 
viscous boundary layer flow will be done using an 


5, — Descriptio n of the calculation method 


integral method formulation of the boundary layer 
equations. With integral methods there is no 
explicit formulation of a turbulence model, but the 
system of equations is supplemented by suitable 
empirical closure relations. The integral method is 
reasonably easy to implement, because in the 
integral method the three-dimensional flow problem 
is formulated as a two-dimensional problem on the 
wing surface and the wake center- surface . This is a 
great advantage in code development. 

Until a few years ago, the boundary layer 
equations have always been used in their steady 
form to compute a steady boundary layer flow. The 
main advantage of this formulation over an unsteady 
formulation has always been its lower computation 
time. However, using the boundary layer equations 
in their unsteady formulation, integrating them 
towards a steady solution, the boundary layer 
method is particularly well suited for 
vectorization and hence for implementation on 
todays supercomputers, see Van Dalsem and Steger 
(Ref. 8), Swafford and Whitfield (Ref. 9). An even 
more important advantage of using the unsteady 
boundary layer equations for solving steady 
boundary layer flow is the robustness of this 
approach. With the unsteady formulation a simple 
time- integration scheme will in any case produce an 
unsteady answer. With the steady formulation a 
space -marching scheme has to be used and the 
specification of initial data proves to be more 
difficult in a space -marching scheme than in a 
time- integration scheme. Also in case the steady 
boundary layer solution is non-unique or non- 
existent, the solution produced by an unsteady 
formulation is probably more useful, than the 
solution produced by a steady formulation. 


Inviscid method 

The full potential equation 
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is solved using a finite volume discretization 

formulated on a curvilinear coordinate system and a 
multigrid method employing ILU/SIP smoothing. At 
present only semi- configurations can be considered. 

The boundary conditions are: 

• on wetted surfaces u„ - qS; (5) 

for the inviscid flow solver, the source 
strength S — 0 on the body, else S is computed 
by the boundary layer solver; 

• in the symmetry plane - 0, (6) 

• in the far-field, except downstream, (p - O’, (7) 

• in the far-field downstream (Trefftz-plane) 
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Two candidate interaction algorithms have been 
considered, namely the semi -Inverse and quasi - 
simultaneous algorithm. (The more sophisticated 
simultaneous algorithm is too expensive to 
Implement in MATRICS because of the fully implicit 
relaxation algorithm employed in the existing 
inviscid potential outer flow solver) . Experience 
by Ashill a. a. (Ref. 10), pp. 35, and Cebeci, Chen 
e.a. (Ref. 11) indicates that the semi -inverse 
algorithm lacks robustness in difficult flow cases. 
Therefore the quasi -simultaneous interaction 
algorithm is used, being the best available 
interaction algorithm that can be implemented in 
interaction with the inviscid potential outer flow 
solution algorithm. 

In the inviscid potential outer flow solver 
the "blowing velocity approach" is used to account 
for the boundary layer effects on the inviscid 
outer flow, which means that a source strength is 
specified at the body surface and a jump in source 
strength at the wake. This approach is reasonably 
standard, while experience by Chen, Li e.a. 
(Ref. 12) reports that this approach is more 
reliable than the boundary layer displacement 
approach . 

Experience by Chow (Ref. 13) indicates that it 
is mandatory to prescribe the pressure jump across 
the wake as a boundary condition to the outer 
inviscid flow. The latter jump Influences the Kutta 
condition prescription in the outer Inviscid flow, 
and this is essential In order to compute realistic 
lift values for the wing. The pressure jump across 
the boundary layer is computed as in Chow (Ref. 13) 
and Lock and Williams (Ref. 14). 


where is the chordwise (wrap-around) grid 
coordinate, 

• across the prescribed vortex sheet 

f tql - q* - q-, (9) 

• (pqS)' - (pqS)- ; 

for the inviscid flow solver the jump across 
the wake, (.)* - (.)", Is computed by the 

boundary layer solver, 

* across the C-0 topology branch cut that 
extends from the tip section to the far-field 
lateral boundary 

[p] * 0, [puj * 0. (10) 


Sul Viscous method 

The steady first order boundary layer 
equations, describing conservation of mass and 
momentum in a general right-handed coordinate 
system, can be found in Myring (Ref. 15). Adding 
time -dependent terms and using first order integral 
thicknesses, the boundary layer equations can be 
integrated with respect to z (normal to the wing 
surface or wake centerline), using the momentum 
equation in z-direction to eliminate the pressure. 
Then the integral equations are obtained in the 
form 


x -momentum: 

1 ££i _ 

q dt dt 


.<r 


r‘-i 


dq 

at 


111 mu i mu ilium i H iminHii i mi 



1 a * u * #„■ 

2-M 2 1 3q .lafJl.u 

h, dx 

hj q 3x J dx \h 1 j 


( 11 ) 


♦ _L£!i2 

h 2 3y 


[lj.au 

[h 2 q dx 




2-M 2 1 3q 

+ iJ.Ci.Kk,. 

h 2 q dy 

J ay [h 2 J 


iii ♦ k 2 i ♦ k,| 

qj % 

q 3y q q 


* ^22^2 



In the latter equation (13) the instationary 
entrainment coefficient is an extension of the 
unsteady two-dimensional definition as used by 
Houwink (Ref. 16). 

Subsequently an expression is needed for the 
calculation of the non-dimensional mass flux S 
through the surface of the wing and wake , 
representing the displacement effect of the 
boundary layer on the inviscid outer flow. Assuming 
an inviscid flow between the stream surface z - 0 
and the displacement surface 5 *, which is a stream 
surface for the inviscid flow, the following 
expression can be obtained from (13) by setting 
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where q is the velocity, J is the Jacobien of the 
transformation from physical to computational 
space, and overbars denote boundary layer edge 
values. The density thickness S fi is the integral of 
(p-p)/p over the boundary layer. 


this Is equivalent to S - w/q, yielding the usual 
Interpretation of S in steady inviscid flow. 

In the used body -conforming non -orthogonal 
coordinate system the y-axis is In the spanwise 
direction of the wing, while the x-axis is in 
chordwise direction wrapping around the wing and 
the wake. 

Next, a streamline coordinate system Is adopted, in 
which the variables (now denoted with tildes) 
reduce to their familiar form, see Myring 
(Ref. 15). Transformation to and from the stream- 
line coordinate system is done whenever necessary. 
Thus It is possible to derive equations using the 
well-known integral parameters in streamline 
coordinates and in the curvilinear system. 

The equations (11) to (14) are thus written in the 
basic variables 7 U , H, q and C, where the cross - 
flow parameter C is defined as 

c „ ai s n (aj)j-«22 (is) 

Reduction of the number of unknowns to the four 
basic variables is established by prescribing a set 
of turbulent velocity profiles in the streamline 
coordinate system, while the density thickness S p 
Is eliminated using the Crocco relation, 
prescribing a parabolic distribution between 
velocity and temperature (Ref. 17). For the time 
being no velocity profiles are used, but proven 
closure relations taken from accepted two- 
dimensional methods for attached as well as 
separated flow (Ref. 16, 18, 19) have been 
implemented in a first version of the code. In a 
next version, more physical closure relations and 
velocity profile families will be used. 

Initial conditions for the turbulent boundary 
layer calculation are generated by the B0LA-2D 
solver (Ref. 20) , which calculates the laminar 
quasi -two -dimensional flow in the leading edge 
region of the wing. 

Boundary conditions are set at the wing root, where 
derivatives in spanwise direction are assumed to be 
zero, and the wing tip, where zero lateral deriva- 
tives in local sweep direction are prescribed, see 
Cross (Ref. 21). Far downstream, a zero gradient 
condition in chordwise direction is specified for 
all quantities. 


The system of equations (11) to (13) is solved 
in combination with an interaction law (to be 
discussed in the next section) . The complete set of 
equations can be written symbolically as 

p 

— Ut + AUj + Buy + Du - f , (16) 
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where u « [0 11( H, q, C] , 

The system (16) appears to be hyperbolic in 
practice. Discretization is done according to the 
directions of the characteristics in (x f t) and 
(y,t) -space, using a matrix-split procedure 
(Ref. 22) . Thus equation (16) is discretized as 


|u, + A* £ x u ♦A'? x u + B* cF y u + B‘? y u + Du * f . 

q 

In smooth parts of the flow second order accurate 
differencing will be obtained using a scheme as for 
example in Spekreijse (Ref. 23). 

The system of equations (17) is solved using the 
fully implicit backward Euler time- integration 
scheme proposed by Steger and Warming (Ref. 24) and 
Yee (Ref. 25). 

U Interaction lav 

In order to avoid a breakdown of the boundary 
layer formulation In separated flow regions an 
extra equation is needed which modifies the 
inviscid flow boundary layer edge velocity q. 
Usually a highly linearized form of the Inviscid 
outer potential flow is taken, for example the two- 
dimensional Hilbert- integral formulation as used by 
Veldman (Ref. 26). An even more simplified form is 
given by Williams (Ref. 27). In this paper the 
latter form is slightly modified, but still derived 
from the linearized potential equation. This will 
be discussed in more detail in section 6.2. In its 
simplest form the interaction law can be written as 



specified on the wing. Further downstream, a quasi - 
simultaneous interaction algorithm is used. In this 
formulation the inviscid flow calculation is done 
in direct mode with a prescribed source strength S 
on the wing and the wake and a prescribed velocity 
jump across the wake. The viscous calculation is 
done in quasi -simultaneous mode with a prescribed 
Inviscid wall velocity, which has been corrected 
for boundary layer curvature effects. Thus the 
boundary layer is computed effectively with a 
prescribed edge velocity instead of the inviscid 
wall velocity, Cebeci, Clark e.a. (Ref. 28) have 
shown that such an approach avoids the initiation 
of undesirable pressure fluctuations in the 
trailing- edge region at reasonably large angles of 
attack. The boundary layer computation computes a 
new source strength S and velocity jump Aq, which 
are used as the subsequent Input for the 

Interactive calculations . 

L. Analysis of Che system of equations 

&U Analysis of boundary layer equation system 

To analyze the properties of the boundary 
layer equations formulated In chapter 5 (Eqs. (1), 
(2), (3), (18)) we assume the following simplifying 
conditions: 

• orthonormal coordinate system, i.e. h 1 -h 2 -l, 

• g "0; k 1 -k 2 -k 3 -0 ; lj— ■ l 2 "l3”0; 

• outer streamline aligned with the x-axis, i.e. 
u/q-i, v/q-0; 

• closure conditions as in Cousteix and 
Houdeville (Ref . 20) , i.e. 

^ 12 " , 8 21 * “ 11 , 

8 22 - -c 2 (« 1 -i 11 ), s 2 - < 19 ) 

With H-Jj/7 u and H 1 -(5-I l )/7 u we obtain an 
equation system 

.5 ♦ Aujj ♦ Buy * f , u « [ ln$ u , H , lnq , InC ] , 

q 

where 


Time -dependent terms are obtained from the 

instationary form of S (equation (4)), paying 
special attention to the limiting case for M 0, 
Two remarks are made: 

1. Equation (18) is written along streamlines, 
which implies that the streamline directions 
are known from the inviscid flow solver and 
are kept fixed during a viscous calculation. 

2. Equation (18) is a law in correction form, 

indicating that it will not influence the 
converged solution. In this way the 

interaction law can be shown to be essential 
to avoid breakdown of the boundary layer 
formulation, but once convergence is obtained 
It does not affect the final solution. 
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5.4 Viscous -inviscid intera ction algorithm 

The leading edge part of the boundary layer 
will be calculated in direct mode using the program 
B0LA-2D (Ref. 20). This presupposes that the flow 
does not separate in this part of the boundary 
layer. The Inviscid flow computation in this part 
Is done In direct mode with a source strength S 
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all values for A are positive for all H, so that 
now all characteristics originate from the upwind 
direction. 

For the case d/dynO we find the following 
characteristic directions in the x,y- plane 
(A-dy/dx) : 

• No interaction law (Ki-O, K^-l) : 


h{ - dHj/dH. 


( 23d > A ■ H-l, 


(27b) 


Setting Krf-1 and Ki-ir/C^Ax) in the third row of 
eqs. (20).. (22) reproduces eq. (18). 

Following Myrlng, reference 16, the 
characteristics of an equation 


H,'A 2 ♦ ((H-UHj-Hid**,) + l+$ r ) A ♦ 

+ 1 ♦ S t - 0; (27c) 


Eu, ♦ Au„ + f - 0, u « [u t . . . uj 
can be obtained from 

det(E-AA) - 0, A - . (24) 

dx 


We first consider the quasi two-dimensional flow 
case (3/3 y— 0) where no interaction law is used, 
i.e. Kj-0, Kj— 1. This way we find the following 
expressions for the characteristic directions A: 


as in the x,t-plane a characteristic direction 
changes from upwind to downwind at separation 
(i.e. at Hjf-0). 

With interaction law (KiH), Kj-1) in case K t -1 , 
M 2 — 0 , €—0 : 


(-(H+1) 2 h(+(H+1)H 1 )A 3 + ( (2H 3 +H z -2H-2)H( ♦ 


♦(-2H 2 *l)H 1 *H 2 +H+l)A 2 ♦ 


A - M,f r (H+l), 


(25b) 


♦ ( (H 2 -1)H,+H 2 +H+2)A +1-0; 


(28b) 


-H(A 2 ♦ ((H+DHj-Hj^l-tH+lJSX) x * 


two positive and one negative value for A are 
found in case of realistic H l (H) -functions . 


- HHj' ♦ H! + (H+l)5 r H 1 / ■ 0. ( 25c ) 

For realistic H^H) - functions , viz. with -H 1( /H 
the roots of the latter quadratic equation are 
real. Consequently, the equation system is fully 
hyperbolic. If H^cO the values for A are greater 
than zero. At separation, that is at H{-0, in the 
minimum of the H^H) curve, one eigenvalue A passes 
through zero, which means that the characteristic 
direction changes from upwind to downwind there. 
Consequently, the equation system models the 
corresponding physical behaviour. 

We continue the analysis for the quasi two- 
dimensional case (3/3 y-0) , but now for the case 
where an Interaction law is used to solve for lnq, 
I.e. K t H0, Ktj-1. In case K t -1 , M 2 — 0 and e-0 we find 
the following expressions for A; 


A - H, 


(26b) 


((H*l) 2 H 1 '-(H*l)H 1 ) A 2 ♦ (-Hj (2H+1) (H+l) ♦ 

♦ H, (2H+3) - 1) A ♦ (H+l) (HH, -Hj) - 0. < 26c > 

For realistic Hj/fi) -functions, namely with (H+l)Hi - 
H^O and HH^ -Hj<0 three values for A are positive, 
one Is negative. For small e<0 It can be shown that 


The characteristics in the x,t-plane for the 
more general case with an Interaction law (K t ^0) 
and without the additional settings K A -1, M^O and 
e-0 can only be analyzed numerically. We find the 
following expressions for the characteristic 
directions : 


A 


H, 


(29a) 


a 3 A 3 ♦ a 2 A 2 + a l A ♦ Sq * 0, (29b) 

where a 3 - det A and Sq - det E. 

A well-defined interaction law has the property 
that a 3 p*0 for all 3. This gives the following 
relation between K t and f^: 

Kt> lUi Kd» 

(HH, -H,) (H+l) (H+l) (1-M 2 ) -H 1 f r (H+l) 

(30) 

if Hi<0 the value K t -0 may be chosen (direct mode). 
If Hi<0 the value ^-0.011^ is suitable for all H 2 
and all ft, in the sense that numerical computation 
of A indicates that then for O'Ol^kl.? all values 
for A are real and positive, while in some rare 
cases complex values are found with a positive real 
part. The boundary layer equation system (and 
interaction law) have consequently favourable 
properties for use in the transonic flow regime at 
cruise conditions. Because all characteristic 
directions generally originate from upstream, 
initial conditions for ? u , H, q and C have to be 



specified upstream, while downstream no initial 
conditions need be specified. 

&- l 2 Properties of Interaction lav 

Following Lock and Williams (Ref. 14), we 
consider a source strength S approximately normal 
to the surface, a velocity q approximately parallel 
to the surface, and define a perturbation potential 
^ as 

B AS -* 8 . (31) 

q 


Parameters are now introduced for the direct and 
inverse parts of eq. (35): 

<3S > 

The parameters K t and ^ can be used to perform 
direct, inverse and quasi -simultaneous computations 
as follows: 

direct : - 0, K* - 1 , (39a) 

inverse ; K t - 1, K* - 0, (39b) 

quasi- simultaneous: K t > 0, K* - 1. (39c) 


The outer potential flow can be described by the 
linearized perturbation potential equation 

(1-M*) - 0, 

with ^ • AS at z ■ 0. (32) 

In the context of the definition of an Interaction 
law we will interprets Aq/q and AS as corrections 
with respect to a starting solution, i.e. 

£2 - 3^1 , AS - S-S*. (33)’ 

q q 

where q k is given by the proceeding outer flow 
computation and S k by the proceeding viscous flow 
calculation. A solution to (32) in subsonic flow in 
Fourier space is. given by 

4 - Ce lM 8 P - ✓jl-M 2 ! . (34) 


With quasi- simultaneous computations, 
inspection of equation (3.8) shows that the 
Inviscld outer flow solution is modified locally at 
places where S differs much from the previous 
solution S*. In the firsts few Iterations between 
inner and outer solver this will lead to non- 
physical velocity distributions in the boundary 
layer, especially at the wake where a velocity 
difference will occur between upper and lower side 
for lifting cases, showing the deficiency of not 
modelling the circulation in the interaction law. 
The inviscld flow solver will therefore have to 
compute this circulation effect on Its own. 


Considering the two-dimensional steady form of 
S and using the property of mass conservation along 
streamlines in the outer flow, we will use In 
eq. (38): 



(40) 


Algebraic manipulation of (31), (32) and (34), 

using the derivative of Aq/q with respect to s in 
order to obtain non- imaginary quantities, then 
leads to a law In Fourier space of the form 

P (iiS-— is!] - ■'(S-S*). (35) 

(q 8s q k 3s J 

In physical space we may consider (35) as the 
leading term In an integral form of an interaction 
law and use It in this simplified form. 

In supersonic flow only the right -running wave is 
considered, 

4 . (36) 

Manipulation now yields a non- imaginary interaction 
law of the form 

0 - - <S-S k ) . (37) 

q 

Considering eq. (37) as the leading term of an 
integral form of an interaction in physical space, 
we now have a form of interaction law that couples 
S to q instead of 3q/3s, which cannot be 
Incorporated In the system of boundary layer 
equations in a simple way. The subsonic law (35) 
appears however to produce useful results both in 
subsonic and supersonic flow, which is mainly due 
to the fact that strong interaction occurs In 
subsonic parts of the flow (at supersonic -subsonic 
shock waves and at trailing edges). This relaxes 
the need for an explicit supersonic lav, also 
considering the findings described in section 5.1. 
We will therefore use eq. (35) both in supersonic 
and subsonic flow. 


It can now be observed that with ^>0 and K^O in 
equation (38) dq/ds has a positive correlation with 
d$*/&s r which is a desired property in subsonic 
flow. 

1. gigllainflrY computational results 

The MATRICS inviscld flow solver has been 
tested extensively and appears to be a robust 
method (Ref. 3, 4). In this section attention will 
be given to the boundary layer solver with 
interaction law and its robustness. Finally a fully 
converged solution of the whole viscous -inviscld 
interaction computation will be shown. 

Starting with an Inviscid velocity distribution 
provided by the MATRICS outer solver, some 
calculations have been made to obtain a converged 
boundary layer solution only. 

Initial values for the boundary layer parameters 
were obtained from a two-dimensional flat plate 
solution. These calculations, without any inter- 
action with the Inviscld outer solver, show the 
changes made by the interaction law. 

The first case is a NACA-0012 straight wing at 
a Mach number of 0.70, zero angle of attack and 
Reynolds number of 9 million. Only the root section 
of the wing will be shown, because this section is 
a symmetry plane with a two-dimensional flow by 
definition, and therefore allows comparison with 
two-dimensional methods. In figure la the inviscld 
velocity distribution at the root section is given. 
It can be shown that the trailing edge stagnation 
behaviour causes a severe boundary layer growth, 
leading to separation if no interaction lav is used 
to adapt the velocity distribution. Figures lb, c, 
d show the resulting velocity, momentum thickness 
and shape factor H at Reynolds number of 9 million. 


calculated with interaction factor Ki-0.05. The 
resulting velocity distribution shows smaller 
decelerating velocity gradients and a less severe 
stagnation behaviour at the trailing edge. In order 
to obtain a physically relevant viscous solution, 
the outflow, computed from this very first boundary 
layer solution, can be added to the full inviscid 
flow solver for an adjustment of the inviscid flow 
to the boundary layer effect. A fully converged 
viscous- inviscid solution requires a number of such 
interations between the inviscid flow solver and 
the boundary layer/interaction law computation. 

The second case is the same straight wing at 
Mach number of 0.799, angle of attack of 2.26 
degrees, and a Reynolds number of 9 million. The 
inviscid velocity distribution at the root section 
(Fig. 2a) will certainly cause separation on the 
upper side of the wing. Starting from a very simple 
flat plate momentum thickness distribution, large 
changes must be expected in the velocity gradients. 
In this case, the resulting momemtum thickness at 
the trailing edge is much below the viscous - 
inviscid converged value, due to the changes in 
velocity distribution, while the shape factor H is 
only showing separation just behind the shock wave 
(Fig. 2b, c, d) . In the viscous- inviscid converged 
solution separation will probably occur from shock 
to trailing edge (Ref. 29) . The non-physical 
velocity jump across the wake (Fig. 2b) is due to 
the asymmetry of the computed flow and the 
simplicity of the interaction law. In the iterative 
process between the full Inviscid flow solver and 
the boundary layer/interaction law, however, the 
errors due to the simplicity of the interaction law 
should disappear, resulting In a converged solution 
and an Inactive interaction law. 

Figures 3a and 3b show the convergence 
histories of the residuals of the x-momentum 
equation, entrainment equation and interaction law 
of the boundary layer equation system for the 
testcases presented in figures 1 and 2. Using the 
same computational parameters, both convergence 
histories show a robust convergence, which means 
that a converged solution can be computed even when 
the inviscid flow velocity distribution and the 
starting solution for the boundary layer flow do 
not at all fit together. 

At the time of writing of the paper, work on 
the Interaction algorithm was making good progress. 
A first fully converged interacting solution has 
been obtained for the attached flow around the 
NACA-0012 wing at a Mach number of 0.70, angle of 
attack of 1.49 degrees, and a Reynolds number of 9 
million. The boundary layer is computed with an 
increasing part of the inviscid flow velocity in 
the first few iterations, while underrelaxation is 
applied to the source strength that is used as 
input to the outer solver. No boundary layer 
curvature effects have yet been accounted for. 
Using this procedure the inner and outer flow are 
smoothly adapted to each other. 

In figure 4 the initial inviscid velocity 
distribution (viz. without a boundary layer effect) 
Is shown, together with the final viscous velocity 
distribution and the corresponding boundary layer 
variables. The difference between the inviscid and 
viscous velocity is small, except at the trailing 
edge, where the boundary layer displacement effect 
is appreciable. This testcase has been calculated 
with thirteen viscous- inviscid iterations in order 
to obtain a converged lift value. The convergence 


history Is shown in figure 5a. Finally, in figure 
5b the convergence history of the residuals In the 
MATRICS outer flow solver are given as a function 
of the number of smoothings. A reduction in both 
maximum and mean value residual of 2.5 orders of 
magnitude has been obtained, which is promising for 
further development of the interaction algorithm. 
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a) Inviscid flow velocity distribution 


b) Boundary layer inviscid edge velocity 



c) Boundary layer streamwiae nonentua thickness d) Boundary layer shape factor H 

Fig. 1 Boundary layer results for root section of a NACA-0012 straight wing at ^*0.70, a*0, Re * 9 * 10® 























a) Inviscid flow velocity distribution 


b) Boundary layer inviscid edge velocity 



c) Boundary layer streasnrise ■oaentua thickness d) Boundary layer shape factor H 


Pig. 2 Boundary layer results for root section of a NACA-0012 straight wing at 1^* 0.799, a* 2.26 
Re » 9 * 10® 





















a) Maximum residual vs. timesteps 
^=0.70, a * 0, Re s 9 * 10 6 


b) Mean residual vs. timesteps 
5^= 0.70, a * 0, Re = 9 *10° 
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c) Maximum residual vs. timesteps 
1^=0.799, a * 2 . 26 , Re * 9 * 10 6 


d) Mean residual vs. timesteps 
1^*0.799, a * 2 . 26 , Re * 9 * 10 e 


Fig. 3 Convergence histories for boundary layer computation on a NACA-0012 straight wing 





c) Boundary layer streamwise ■omentum thickness d) Boundary layer shape factor H 


Fig. 4 Viscous-inviscid interaction results for root section of a NACA-0012 straight wing at ^*0.70, 
a » 1 . 49 , Re * 9 * 10 6 
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